home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / pcl / src-16f.lha / code / float.lisp < prev    next >
Encoding:
Text File  |  1992-05-30  |  29.4 KB  |  845 lines
  1. ;;; -*- Mode: Lisp; Package: KERNEL; Log: code.log -*-
  2. ;;;
  3. ;;; **********************************************************************
  4. ;;; This code was written as part of the CMU Common Lisp project at
  5. ;;; Carnegie Mellon University, and has been placed in the public domain.
  6. ;;; If you want to use this code or any part of CMU Common Lisp, please contact
  7. ;;; Scott Fahlman or slisp-group@cs.cmu.edu.
  8. ;;;
  9. (ext:file-comment
  10.   "$Header: float.lisp,v 1.9 92/02/07 11:20:16 ram Exp $")
  11. ;;;
  12. ;;; **********************************************************************
  13. ;;;
  14. ;;;    This file contains the definitions of float specific number support
  15. ;;; (other than irrational stuff, which is in irrat.)  There is code in here
  16. ;;; that assumes there are only two float formats: IEEE single and double.
  17. ;;;
  18. ;;; Author: Rob MacLachlan
  19. ;;; 
  20. (in-package "KERNEL")
  21. (export '(%unary-truncate %unary-round))
  22.  
  23. (in-package "LISP")
  24. (export '(least-positive-normalized-short-float
  25.       least-positive-normalized-single-float
  26.       least-positive-normalized-double-float
  27.       least-positive-normalized-long-float
  28.       least-negative-normalized-short-float
  29.       least-negative-normalized-single-float
  30.       least-negative-normalized-double-float
  31.       least-negative-normalized-long-float
  32.       least-positive-single-float
  33.       least-positive-short-float
  34.       least-negative-single-float
  35.       least-negative-short-float
  36.       least-positive-double-float
  37.       least-positive-long-float
  38.       least-negative-double-float
  39.       least-negative-long-float
  40.       most-positive-single-float
  41.       most-positive-short-float
  42.       most-negative-single-float
  43.       most-negative-short-float
  44.       most-positive-double-float
  45.       most-positive-long-float
  46.       most-negative-double-float
  47.       most-negative-long-float))
  48.  
  49. (in-package "EXTENSIONS")
  50. (export '(single-float-positive-infinity short-float-positive-infinity
  51.       double-float-positive-infinity long-float-positive-infinity
  52.       single-float-negative-infinity short-float-negative-infinity
  53.       double-float-negative-infinity long-float-negative-infinity
  54.       set-floating-point-modes float-denormalized-p float-nan-p
  55.       float-trapping-nan-p float-infinity-p))
  56.  
  57. (in-package "KERNEL")
  58.  
  59.  
  60. ;;;; Utilities:
  61.  
  62. ;;; SINGLE-FROM-BITS, DOUBLE-FROM-BITS  --  Internal
  63. ;;;
  64. ;;;    These functions let us create floats from bits with the significand
  65. ;;; uniformly represented as an integer.  This is less efficient for double
  66. ;;; floats, but is more convenient when making special values, etc.
  67. ;;;
  68. (defun single-from-bits (sign exp sig)
  69.   (declare (type bit sign) (type (unsigned-byte 24) sig)
  70.        (type (unsigned-byte 8) exp))
  71.   (make-single-float
  72.    (dpb exp vm:single-float-exponent-byte
  73.     (dpb sig vm:single-float-significand-byte
  74.          (if (zerop sign) 0 -1)))))
  75. ;;;
  76. (defun double-from-bits (sign exp sig)
  77.   (declare (type bit sign) (type (unsigned-byte 53) sig)
  78.        (type (unsigned-byte 11) exp))
  79.   (make-double-float (dpb exp vm:double-float-exponent-byte
  80.               (dpb (ash sig -32) vm:double-float-significand-byte
  81.                    (if (zerop sign) 0 -1)))
  82.              (ldb (byte 32 0) sig)))
  83.                     
  84.  
  85. ;;;; Float parameters:
  86.  
  87. (defconstant least-positive-single-float (single-from-bits 0 0 1))
  88. (defconstant least-positive-short-float least-positive-single-float)
  89. (defconstant least-negative-single-float (single-from-bits 1 0 1))
  90. (defconstant least-negative-short-float least-negative-single-float)
  91. (defconstant least-positive-double-float (double-from-bits 0 0 1))
  92. (defconstant least-positive-long-float least-positive-double-float)
  93. (defconstant least-negative-double-float (double-from-bits 1 0 1))
  94. (defconstant least-negative-long-float least-negative-double-float)
  95.  
  96. (defconstant least-positive-normalized-single-float
  97.   (single-from-bits 0 vm:single-float-normal-exponent-min 0))
  98. (defconstant least-positive-normalized-short-float
  99.   least-positive-normalized-single-float)
  100. (defconstant least-negative-normalized-single-float
  101.   (single-from-bits 1 vm:single-float-normal-exponent-min 0))
  102. (defconstant least-negative-normalized-short-float
  103.   least-negative-normalized-single-float)
  104. (defconstant least-positive-normalized-double-float
  105.   (double-from-bits 0 vm:double-float-normal-exponent-min 0))
  106. (defconstant least-positive-normalized-long-float
  107.   least-positive-normalized-double-float)
  108. (defconstant least-negative-normalized-double-float
  109.   (double-from-bits 1 vm:double-float-normal-exponent-min 0))
  110. (defconstant least-negative-normalized-long-float
  111.   least-negative-normalized-double-float)
  112.  
  113. (defconstant most-positive-single-float
  114.   (single-from-bits 0 vm:single-float-normal-exponent-max
  115.             (ldb vm:single-float-significand-byte -1)))
  116. (defconstant most-positive-short-float most-positive-single-float)
  117. (defconstant most-negative-single-float
  118.   (single-from-bits 1 vm:single-float-normal-exponent-max
  119.             (ldb vm:single-float-significand-byte -1)))
  120. (defconstant most-negative-short-float most-negative-single-float)
  121. (defconstant most-positive-double-float
  122.   (double-from-bits 0 vm:double-float-normal-exponent-max
  123.             (ldb vm:double-float-significand-byte -1)))
  124. (defconstant most-positive-long-float most-positive-double-float)
  125. (defconstant most-negative-double-float
  126.   (double-from-bits 1 vm:double-float-normal-exponent-max
  127.             (ldb vm:double-float-significand-byte -1)))
  128. (defconstant most-negative-long-float most-negative-double-float)
  129.  
  130. (defconstant single-float-positive-infinity
  131.   (single-from-bits 0 (1+ vm:single-float-normal-exponent-max) 0))
  132. (defconstant short-float-positive-infinity single-float-positive-infinity)
  133. (defconstant single-float-negative-infinity
  134.   (single-from-bits 1 (1+ vm:single-float-normal-exponent-max) 0))
  135. (defconstant short-float-negative-infinity single-float-negative-infinity)
  136. (defconstant double-float-positive-infinity
  137.   (double-from-bits 0 (1+ vm:double-float-normal-exponent-max) 0))
  138. (defconstant long-float-positive-infinity double-float-positive-infinity)
  139. (defconstant double-float-negative-infinity
  140.   (double-from-bits 1 (1+ vm:double-float-normal-exponent-max) 0))
  141. (defconstant long-float-negative-infinity double-float-negative-infinity)
  142.  
  143. (defconstant single-float-epsilon
  144.   (single-from-bits 0 (- vm:single-float-bias (1- vm:single-float-digits)) 1))
  145. (defconstant short-float-epsilon single-float-epsilon)
  146. (defconstant single-float-negative-epsilon
  147.   (single-from-bits 0 (- vm:single-float-bias vm:single-float-digits) 1))
  148. (defconstant short-float-negative-epsilon single-float-negative-epsilon)
  149. (defconstant double-float-epsilon
  150.   (double-from-bits 0 (- vm:double-float-bias (1- vm:double-float-digits)) 1))
  151. (defconstant long-float-epsilon double-float-epsilon)
  152. (defconstant double-float-negative-epsilon
  153.   (double-from-bits 0 (- vm:double-float-bias vm:double-float-digits) 1))
  154. (defconstant long-float-negative-epsilon double-float-negative-epsilon)
  155.  
  156.  
  157. ;;;; Float predicates and environment query:
  158.  
  159. (proclaim '(maybe-inline float-denormalized-p float-infinity-p float-nan-p
  160.              float-trapping-nan-p))
  161.  
  162. ;;; FLOAT-DENORMALIZED-P  --  Public
  163. ;;;
  164. (defun float-denormalized-p (x)
  165.   "Return true if the float X is denormalized."
  166.   (number-dispatch ((x float))
  167.     ((single-float)
  168.      (and (zerop (ldb vm:single-float-exponent-byte (single-float-bits x)))
  169.       (not (zerop x))))
  170.     ((double-float)
  171.      (and (zerop (ldb vm:double-float-exponent-byte
  172.               (double-float-high-bits x)))
  173.       (not (zerop x))))))
  174.  
  175. (macrolet ((frob (name doc single double)
  176.          `(defun ,name (x)
  177.         ,doc
  178.         (number-dispatch ((x float))
  179.           ((single-float)
  180.            (let ((bits (single-float-bits x)))
  181.              (and (> (ldb vm:single-float-exponent-byte bits)
  182.                  vm:single-float-normal-exponent-max)
  183.               ,single)))
  184.           ((double-float)
  185.            (let ((hi (double-float-high-bits x))
  186.              (lo (double-float-low-bits x)))
  187.              (and (> (ldb vm:double-float-exponent-byte hi)
  188.                  vm:double-float-normal-exponent-max)
  189.               ,double)))))))
  190.  
  191.   (frob float-infinity-p "Return true if the float X is an infinity (+ or -)."
  192.     (zerop (ldb vm:single-float-significand-byte bits))
  193.     (and (zerop (ldb vm:double-float-significand-byte hi))
  194.      (zerop lo)))
  195.  
  196.   (frob float-nan-p "Return true if the float X is a NaN (Not a Number)."
  197.     (not (zerop (ldb vm:single-float-significand-byte bits)))
  198.     (or (not (zerop (ldb vm:double-float-significand-byte hi)))
  199.     (not (zerop lo))))
  200.  
  201.   (frob float-trapping-nan-p
  202.     "Return true if the float X is a trapping NaN (Not a Number)."
  203.     (not (zerop (logand (ldb vm:single-float-significand-byte bits)
  204.             vm:single-float-trapping-nan-bit)))
  205.     (progn
  206.       lo; ignore
  207.       (not (zerop (logand (ldb vm:double-float-significand-byte hi)
  208.               vm:double-float-trapping-nan-bit))))))
  209.  
  210.  
  211. ;;; FLOAT-PRECISION  --  Public
  212. ;;;
  213. ;;;    If denormalized, use a subfunction from INTEGER-DECODE-FLOAT to find the
  214. ;;; actual exponent (and hence how denormalized it is), otherwise we just
  215. ;;; return the number of digits or 0.
  216. ;;;
  217. (proclaim '(maybe-inline float-precision))
  218. (defun float-precision (f)
  219.   "Returns a non-negative number of significant digits in it's float argument.
  220.   Will be less than FLOAT-DIGITS if denormalized or zero."
  221.   (macrolet ((frob (digits bias decode)
  222.            `(cond ((zerop f) 0)
  223.               ((float-denormalized-p f)
  224.                (multiple-value-bind (ignore exp)
  225.                         (,decode f)
  226.              (declare (ignore ignore))
  227.              (truly-the fixnum
  228.                     (+ ,digits (1- ,digits) ,bias exp))))
  229.               (t
  230.                ,digits))))
  231.     (number-dispatch ((f float))
  232.       ((single-float)
  233.        (frob vm:single-float-digits vm:single-float-bias
  234.      integer-decode-single-denorm))
  235.       ((double-float)
  236.        (frob vm:double-float-digits vm:double-float-bias
  237.      integer-decode-double-denorm)))))
  238.  
  239.  
  240. (defun float-sign (float1 &optional (float2 (float 1 float1)))
  241.   "Returns a floating-point number that has the same sign as
  242.    float1 and, if float2 is given, has the same absolute value
  243.    as float2."
  244.   (declare (float float1 float2))
  245.   (float-sign float1 float2))
  246.  
  247. (defun float-format-digits (format)
  248.   (ecase format
  249.     ((short-float single-float) vm:single-float-digits)
  250.     ((double-float long-float) vm:double-float-digits)))
  251.  
  252. (proclaim '(inline float-digits float-radix))
  253.  
  254. (defun float-digits (f)
  255.   "Returns a non-negative number of radix-b digits used in the
  256.    representation of it's argument.  See Common Lisp: The Language
  257.    by Guy Steele for more details."
  258.   (number-dispatch ((f float))
  259.     ((single-float) vm:single-float-digits)
  260.     ((double-float) vm:double-float-digits)))
  261.  
  262. (defun float-radix (f)
  263.   "Returns (as an integer) the radix b of its floating-point
  264.    argument."
  265.   (declare (ignore f))
  266.   2)
  267.  
  268.  
  269.  
  270. ;;;; INTEGER-DECODE-FLOAT and DECODE-FLOAT:
  271.  
  272. (proclaim '(maybe-inline integer-decode-single-float
  273.              integer-decode-double-float))
  274.  
  275. ;;; INTEGER-DECODE-SINGLE-DENORM  --  Internal
  276. ;;;
  277. ;;;    Handle the denormalized case of INTEGER-DECODE-FLOAT for SINGLE-FLOAT.
  278. ;;;
  279. (defun integer-decode-single-denorm (x)
  280.   (declare (type single-float x))
  281.   (let* ((bits (single-float-bits (abs x)))
  282.      (sig (ash (ldb vm:single-float-significand-byte bits) 1))
  283.      (extra-bias 0))
  284.     (declare (type (unsigned-byte 24) sig)
  285.          (type (integer 0 23) extra-bias))
  286.     (loop
  287.       (unless (zerop (logand sig vm:single-float-hidden-bit))
  288.     (return))
  289.       (setq sig (ash sig 1))
  290.       (incf extra-bias))
  291.     (values sig
  292.         (- (- vm:single-float-bias) vm:single-float-digits extra-bias)
  293.         (if (minusp (float-sign x)) -1 1))))
  294.  
  295.  
  296. ;;; INTEGER-DECODE-SINGLE-FLOAT  --  Internal
  297. ;;;
  298. ;;;    Handle the single-float case of INTEGER-DECODE-FLOAT.  If an infinity or
  299. ;;; NAN, error.  If a denorm, call i-d-s-DENORM to handle it.
  300. ;;;
  301. (defun integer-decode-single-float (x)
  302.   (declare (single-float x))
  303.   (let* ((bits (single-float-bits (abs x)))
  304.      (exp (ldb vm:single-float-exponent-byte bits))
  305.      (sig (ldb vm:single-float-significand-byte bits))
  306.      (sign (if (minusp (float-sign x)) -1 1))
  307.      (biased (- exp vm:single-float-bias vm:single-float-digits)))
  308.     (declare (fixnum biased))
  309.     (unless (<= exp vm:single-float-normal-exponent-max)
  310.       (error "Can't decode NAN or infinity: ~S." x))
  311.     (cond ((and (zerop exp) (zerop sig))
  312.        (values 0 biased sign))
  313.       ((< exp vm:single-float-normal-exponent-min)
  314.        (integer-decode-single-denorm x))
  315.       (t
  316.        (values (logior sig vm:single-float-hidden-bit) biased sign)))))
  317.  
  318.  
  319. ;;; INTEGER-DECODE-DOUBLE-DENORM  --  Internal
  320. ;;;
  321. ;;;    Like INTEGER-DECODE-SINGLE-DENORM, only doubly so.
  322. ;;;
  323. (defun integer-decode-double-denorm (x)
  324.   (declare (type double-float x))
  325.   (let* ((high-bits (double-float-high-bits (abs x)))
  326.      (sig-high (ldb vm:double-float-significand-byte high-bits))
  327.      (low-bits (double-float-low-bits x))
  328.      (sign (if (minusp (float-sign x)) -1 1))
  329.      (biased (- (- vm:double-float-bias) vm:double-float-digits)))
  330.     (if (zerop sig-high)
  331.     (let ((sig low-bits)
  332.           (extra-bias (- vm:double-float-digits 33))
  333.           (bit (ash 1 31)))
  334.       (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
  335.       (loop
  336.         (unless (zerop (logand sig bit)) (return))
  337.         (setq sig (ash sig 1))
  338.         (incf extra-bias))
  339.       (values (ash sig (- vm:double-float-digits 32))
  340.           (truly-the fixnum (- biased extra-bias))
  341.           sign))
  342.     (let ((sig (ash sig-high 1))
  343.           (extra-bias 0))
  344.       (declare (type (unsigned-byte 32) sig) (fixnum extra-bias))
  345.       (loop
  346.         (unless (zerop (logand sig vm:double-float-hidden-bit))
  347.           (return))
  348.         (setq sig (ash sig 1))
  349.         (incf extra-bias))
  350.       (values (logior (ash sig 32) (ash low-bits (1- extra-bias)))
  351.           (truly-the fixnum (- biased extra-bias))
  352.           sign)))))
  353.  
  354.  
  355. ;;; INTEGER-DECODE-DOUBLE-FLOAT  --  Internal
  356. ;;;
  357. ;;;    Like INTEGER-DECODE-SINGLE-FLOAT, only doubly so.
  358. ;;;
  359. (defun integer-decode-double-float (x)
  360.   (declare (double-float x))
  361.   (let* ((abs (abs x))
  362.      (hi (double-float-high-bits abs))
  363.      (lo (double-float-low-bits abs))
  364.      (exp (ldb vm:double-float-exponent-byte hi))
  365.      (sig (ldb vm:double-float-significand-byte hi))
  366.      (sign (if (minusp (float-sign x)) -1 1))
  367.      (biased (- exp vm:double-float-bias vm:double-float-digits)))
  368.     (declare (fixnum biased))
  369.     (unless (<= exp vm:double-float-normal-exponent-max)
  370.       (error "Can't decode NAN or infinity: ~S." x))
  371.     (cond ((and (zerop exp) (zerop sig) (zerop lo))
  372.        (values 0 biased sign))
  373.       ((< exp vm:double-float-normal-exponent-min)
  374.        (integer-decode-double-denorm x))
  375.       (t
  376.        (values
  377.         (logior (ash (logior (ldb vm:double-float-significand-byte hi)
  378.                  vm:double-float-hidden-bit)
  379.              32)
  380.             lo)
  381.         biased sign)))))
  382.  
  383.  
  384. ;;; INTEGER-DECODE-FLOAT  --  Public
  385. ;;;
  386. ;;;    Dispatch to the correct type-specific i-d-f function.
  387. ;;;
  388. (defun integer-decode-float (x)
  389.   "Returns three values:
  390.    1) an integer representation of the significand.
  391.    2) the exponent for the power of 2 that the significand must be multiplied
  392.       by to get the actual value.  This differs from the DECODE-FLOAT exponent
  393.       by FLOAT-DIGITS, since the significand has been scaled to have all its
  394.       digits before the radix point.
  395.    3) -1 or 1 (i.e. the sign of the argument.)"
  396.   (number-dispatch ((x float))
  397.     ((single-float)
  398.      (integer-decode-single-float x))
  399.     ((double-float)
  400.      (integer-decode-double-float x))))
  401.  
  402.  
  403. (proclaim '(maybe-inline decode-single-float decode-double-float))
  404.  
  405. ;;; DECODE-SINGLE-DENORM  --  Internal
  406. ;;;
  407. ;;;    Handle the denormalized case of DECODE-SINGLE-FLOAT.  We call
  408. ;;; INTEGER-DECODE-SINGLE-DENORM and then make the result into a float.
  409. ;;;
  410. (defun decode-single-denorm (x)
  411.   (declare (type single-float x))
  412.   (multiple-value-bind (sig exp sign)
  413.                (integer-decode-single-denorm x)
  414.     (values (make-single-float
  415.          (dpb sig vm:single-float-significand-byte
  416.           (dpb vm:single-float-bias vm:single-float-exponent-byte 0)))
  417.         (truly-the fixnum (+ exp vm:single-float-digits))
  418.         (float sign x))))
  419.  
  420.  
  421. ;;; DECODE-SINGLE-FLOAT  --  Internal
  422. ;;;
  423. ;;;    Handle the single-float case of DECODE-FLOAT.  If an infinity or NAN,
  424. ;;; error.  If a denorm, call d-s-DENORM to handle it.
  425. ;;;
  426. (defun decode-single-float (x)
  427.   (declare (single-float x))
  428.   (let* ((bits (single-float-bits (abs x)))
  429.      (exp (ldb vm:single-float-exponent-byte bits))
  430.      (sign (float-sign x))
  431.      (biased (truly-the single-float-exponent
  432.                 (- exp vm:single-float-bias))))
  433.     (unless (<= exp vm:single-float-normal-exponent-max) 
  434.       (error "Can't decode NAN or infinity: ~S." x))
  435.     (cond ((zerop x)
  436.        (values 0.0f0 biased sign))
  437.       ((< exp vm:single-float-normal-exponent-min)
  438.        (decode-single-denorm x))
  439.       (t
  440.        (values (make-single-float
  441.             (dpb vm:single-float-bias
  442.              vm:single-float-exponent-byte
  443.              bits))
  444.            biased sign)))))
  445.  
  446.  
  447. ;;; DECODE-DOUBLE-DENORM  --  Internal
  448. ;;;
  449. ;;;    Like DECODE-SINGLE-DENORM, only doubly so.
  450. ;;; 
  451. (defun decode-double-denorm (x)
  452.   (declare (double-float x))
  453.   (multiple-value-bind (sig exp sign)
  454.                (integer-decode-double-denorm x)
  455.     (values (make-double-float
  456.          (dpb (logand (ash sig -32) (lognot vm:double-float-hidden-bit))
  457.           vm:double-float-significand-byte
  458.           (dpb vm:double-float-bias vm:double-float-exponent-byte 0))
  459.          (ldb (byte 32 0) sig))
  460.         (truly-the fixnum (+ exp vm:double-float-digits))
  461.         (float sign x))))
  462.  
  463.  
  464. ;;; DECODE-DOUBLE-FLOAT  --  Public
  465. ;;;
  466. ;;;    Like DECODE-SINGLE-FLOAT, only doubly so.
  467. ;;;
  468. (defun decode-double-float (x)
  469.   (declare (double-float x))
  470.   (let* ((abs (abs x))
  471.      (hi (double-float-high-bits abs))
  472.      (lo (double-float-low-bits abs))
  473.      (exp (ldb vm:double-float-exponent-byte hi))
  474.      (sign (float-sign x))
  475.      (biased (truly-the double-float-exponent
  476.                 (- exp vm:double-float-bias))))
  477.     (unless (<= exp vm:double-float-normal-exponent-max)
  478.       (error "Can't decode NAN or infinity: ~S." x))
  479.     (cond ((zerop x)
  480.        (values 0.0d0 biased sign))
  481.       ((< exp vm:double-float-normal-exponent-min)
  482.        (decode-double-denorm x))
  483.       (t
  484.        (values (make-double-float
  485.             (dpb vm:double-float-bias vm:double-float-exponent-byte hi)
  486.             lo)
  487.            biased sign)))))
  488.  
  489.  
  490. ;;; DECODE-FLOAT  --  Public
  491. ;;;
  492. ;;;    Dispatch to the appropriate type-specific function.
  493. ;;;
  494. (defun decode-float (f)
  495.   "Returns three values:
  496.    1) a floating-point number representing the significand.  This is always
  497.       between 0.5 (inclusive) and 1.0 (exclusive).
  498.    2) an integer representing the exponent.
  499.    3) -1.0 or 1.0 (i.e. the sign of the argument.)"
  500.   (number-dispatch ((f float))
  501.     ((single-float)
  502.      (decode-single-float f))
  503.     ((double-float)
  504.      (decode-double-float f))))
  505.  
  506.  
  507. ;;;; SCALE-FLOAT:
  508.  
  509. (proclaim '(maybe-inline scale-single-float scale-double-float))
  510.  
  511. ;;; SCALE-FLOAT-MAYBE-UNDERFLOW  --  Internal
  512. ;;;
  513. ;;;    Handle float scaling where the X is denormalized or the result is
  514. ;;; denormalized or underflows to 0.
  515. ;;;
  516. (defun scale-float-maybe-underflow (x exp)
  517.   (multiple-value-bind (sig old-exp)
  518.                (integer-decode-float x)
  519.     (let* ((digits (float-digits x))
  520.        (new-exp (+ exp old-exp digits
  521.                (etypecase x
  522.              (single-float vm:single-float-bias)
  523.              (double-float vm:double-float-bias))))
  524.        (sign (if (minusp (float-sign x)) 1 0)))
  525.       (cond
  526.        ((< new-exp
  527.        (etypecase x
  528.          (single-float vm:single-float-normal-exponent-min)
  529.          (double-float vm:double-float-normal-exponent-min)))
  530.     (when (vm:current-float-trap :inexact)
  531.       (error 'floating-point-inexact :operation 'scale-float
  532.          :operands (list x exp)))
  533.     (when (vm:current-float-trap :underflow)
  534.       (error 'floating-point-underflow :operation 'scale-float
  535.          :operands (list x exp)))
  536.     (let ((shift (1- new-exp)))
  537.       (if (< shift (- (1- digits)))
  538.           (float-sign x 0.0)
  539.           (etypecase x
  540.         (single-float (single-from-bits sign 0 (ash sig shift)))
  541.         (double-float (double-from-bits sign 0 (ash sig shift)))))))
  542.        (t
  543.     (etypecase x
  544.       (single-float (single-from-bits sign new-exp sig))
  545.       (double-float (double-from-bits sign new-exp sig))))))))
  546.  
  547.  
  548. ;;; SCALE-FLOAT-MAYBE-OVERFLOW  --  Internal
  549. ;;;
  550. ;;;    Called when scaling a float overflows, or the oringinal float was a NaN
  551. ;;; or infinity.  If overflow errors are trapped, then error, otherwise return
  552. ;;; the appropriate infinity.  If a NaN, signal or not as appropriate.
  553. ;;;
  554. (defun scale-float-maybe-overflow (x exp)
  555.   (cond
  556.    ((float-infinity-p x)
  557.     ;; Infinity is infinity, no matter how small...
  558.     x)
  559.    ((float-nan-p x)
  560.     (when (and (float-trapping-nan-p x)
  561.            (vm:current-float-trap :invalid))
  562.       (error 'floating-point-invalid :operation 'scale-float
  563.          :operands (list x exp)))
  564.     x)
  565.    (t
  566.     (when (vm:current-float-trap :overflow)
  567.       (error 'floating-point-overflow :operation 'scale-float
  568.          :operands (list x exp)))
  569.     (when (vm:current-float-trap :inexact)
  570.       (error 'floating-point-inexact :operation 'scale-float
  571.          :operands (list x exp)))
  572.     (* (float-sign x)
  573.        (etypecase x
  574.      (single-float single-float-positive-infinity)
  575.      (double-float double-float-positive-infinity))))))
  576.  
  577.  
  578. ;;; SCALE-SINGLE-FLOAT, SCALE-DOUBLE-FLOAT  --  Internal
  579. ;;;
  580. ;;;    Scale a single or double float, calling the correct over/underflow
  581. ;;; functions.
  582. ;;;
  583. (defun scale-single-float (x exp)
  584.   (declare (single-float x) (fixnum exp))
  585.   (let* ((bits (single-float-bits x))
  586.      (old-exp (ldb vm:single-float-exponent-byte bits))
  587.      (new-exp (+ old-exp exp)))
  588.     (cond
  589.      ((zerop x) x)
  590.      ((or (< old-exp vm:single-float-normal-exponent-min)
  591.       (< new-exp vm:single-float-normal-exponent-min))
  592.       (scale-float-maybe-underflow x exp))
  593.      ((or (> old-exp vm:single-float-normal-exponent-max)
  594.       (> new-exp vm:single-float-normal-exponent-max))
  595.       (scale-float-maybe-overflow x exp))
  596.      (t
  597.       (make-single-float (dpb new-exp vm:single-float-exponent-byte bits))))))
  598. ;;;
  599. (defun scale-double-float (x exp)
  600.   (declare (double-float x) (fixnum exp))
  601.   (let* ((hi (double-float-high-bits x))
  602.      (lo (double-float-low-bits x))
  603.      (old-exp (ldb vm:double-float-exponent-byte hi))
  604.      (new-exp (+ old-exp exp)))
  605.     (cond
  606.      ((zerop x) x)
  607.      ((or (< old-exp vm:double-float-normal-exponent-min)
  608.       (< new-exp vm:double-float-normal-exponent-min))
  609.       (scale-float-maybe-underflow x exp))
  610.      ((or (> old-exp vm:double-float-normal-exponent-max)
  611.       (> new-exp vm:double-float-normal-exponent-max))
  612.       (scale-float-maybe-overflow x exp))
  613.      (t
  614.       (make-double-float (dpb new-exp vm:double-float-exponent-byte hi)
  615.              lo)))))
  616.  
  617.  
  618. ;;; SCALE-FLOAT  --  Public
  619. ;;;
  620. ;;;    Dispatch to the correct type-specific scale-float function.
  621. ;;;
  622. (defun scale-float (f ex)
  623.   "Returns the value (* f (expt (float 2 f) ex)), but with no unnecessary loss
  624.   of precision or overflow."
  625.   (number-dispatch ((f float))
  626.     ((single-float)
  627.      (scale-single-float f ex))
  628.     ((double-float)
  629.      (scale-double-float f ex))))
  630.  
  631.  
  632. ;;;; Converting to/from floats:
  633.  
  634. (defun float (number &optional (other () otherp))
  635.   "Converts any REAL to a float.  If OTHER is not provided, it returns a
  636.   SINGLE-FLOAT if NUMBER is not already a FLOAT.  If OTHER is provided, the
  637.   result is the same float format as OTHER."
  638.   (if otherp
  639.       (number-dispatch ((number real) (other float))
  640.     (((foreach rational single-float double-float)
  641.       (foreach single-float double-float))
  642.      (coerce number '(dispatch-type other))))
  643.       (if (floatp number)
  644.       number
  645.       (coerce number 'single-float))))
  646.  
  647.  
  648. (macrolet ((frob (name type)
  649.          `(defun ,name (x)
  650.         (number-dispatch ((x real))
  651.           (((foreach single-float double-float fixnum))
  652.            (coerce x ',type))
  653.           ((bignum)
  654.            (bignum-to-float x ',type))
  655.           ((ratio)
  656.            (let ((num (numerator x))
  657.              (den (denominator x)))
  658.              (if (and (fixnump num) (fixnump den))
  659.              (/ (coerce num ',type) (coerce den ',type))
  660.              (float-bignum-ratio x ',type))))))))
  661.   (frob %single-float single-float)
  662.   (frob %double-float double-float))
  663.  
  664.  
  665. #|
  666. These might be useful if we ever have a machine w/o float/integer conversion
  667. hardware.  For now, we'll use special ops that uninterruptibly frob the
  668. rounding modes & do ieee round-to-integer.
  669.  
  670. ;;; %UNARY-TRUNCATE-SINGLE-FLOAT/FIXNUM  --  Interface
  671. ;;;
  672. ;;;    The compiler compiles a call to this when we are doing %UNARY-TRUNCATE
  673. ;;; and the result is known to be a fixnum.  We can avoid some generic
  674. ;;; arithmetic in this case.
  675. ;;;
  676. (defun %unary-truncate-single-float/fixnum (x)
  677.   (declare (single-float x) (values fixnum))
  678.   (locally (declare (optimize (speed 3) (safety 0)))
  679.     (let* ((bits (single-float-bits x))
  680.        (exp (ldb vm:single-float-exponent-byte bits))
  681.        (frac (logior (ldb vm:single-float-significand-byte bits)
  682.              vm:single-float-hidden-bit))
  683.        (shift (- exp vm:single-float-digits vm:single-float-bias)))
  684.       (when (> exp vm:single-float-normal-exponent-max)
  685.     (error 'floating-point-invalid :operator 'truncate
  686.            :operands (list x)))
  687.       (if (<= shift (- vm:single-float-digits))
  688.       0
  689.       (let ((res (ash frac shift)))
  690.         (declare (type (unsigned-byte 31) res)) 
  691.         (if (minusp bits)
  692.         (- res)
  693.         res))))))
  694.  
  695.  
  696. ;;; %UNARY-TRUNCATE-DOUBLE-FLOAT/FIXNUM  --  Interface
  697. ;;;
  698. ;;;    Double-float version of this operation (see above single op).
  699. ;;;
  700. (defun %unary-truncate-double-float/fixnum (x)
  701.   (declare (double-float x) (values fixnum))
  702.   (locally (declare (optimize (speed 3) (safety 0)))
  703.     (let* ((hi-bits (double-float-high-bits x))
  704.        (exp (ldb vm:double-float-exponent-byte hi-bits))
  705.        (frac (logior (ldb vm:double-float-significand-byte hi-bits)
  706.              vm:double-float-hidden-bit))
  707.        (shift (- exp (- vm:double-float-digits vm:word-bits)
  708.              vm:double-float-bias)))
  709.       (when (> exp vm:double-float-normal-exponent-max)
  710.     (error 'floating-point-invalid :operator 'truncate
  711.            :operands (list x)))
  712.       (if (<= shift (- vm:word-bits vm:double-float-digits))
  713.       0
  714.       (let* ((res-hi (ash frac shift))
  715.          (res (if (plusp shift)
  716.               (logior res-hi
  717.                   (the fixnum
  718.                        (ash (double-float-low-bits x)
  719.                         (- shift vm:word-bits))))
  720.               res-hi)))
  721.         (declare (type (unsigned-byte 31) res-hi res))
  722.         (if (minusp hi-bits)
  723.         (- res)
  724.         res))))))
  725. |#
  726.  
  727.   
  728. ;;; %UNARY-TRUNCATE  --  Interface
  729. ;;;
  730. ;;;    This function is called when we are doing a truncate without any funky
  731. ;;; divisor, i.e. converting a float or ratio to an integer.  Note that we do
  732. ;;; *not* return the second value of truncate, so it must be computed by the
  733. ;;; caller if needed.
  734. ;;;
  735. ;;;    In the float case, we pick off small arguments so that compiler can use
  736. ;;; special-case operations.  We use an exclusive test, since (due to round-off
  737. ;;; error), (float most-positive-fixnum) may be greater than
  738. ;;; most-positive-fixnum.
  739. ;;;
  740. (defun %unary-truncate (number)
  741.   (number-dispatch ((number real))
  742.     ((integer) number)
  743.     ((ratio) (values (truncate (numerator number) (denominator number))))
  744.     (((foreach single-float double-float))
  745.      (if (< (float most-negative-fixnum number)
  746.         number
  747.         (float most-positive-fixnum number))
  748.      (truly-the fixnum (%unary-truncate number))
  749.      (multiple-value-bind (bits exp)
  750.                   (integer-decode-float number)
  751.        (let ((res (ash bits exp)))
  752.          (if (minusp number)
  753.          (- res)
  754.          res)))))))
  755.  
  756.  
  757. ;;; %UNARY-ROUND  --  Interface
  758. ;;;
  759. ;;;    Similar to %UNARY-TRUNCATE, but rounds to the nearest integer.  If we
  760. ;;; can't use the round primitive, then we do our own round-to-nearest on the
  761. ;;; result of i-d-f.  [Note that this rounding will really only happen with
  762. ;;; double floats, since the whole single-float fraction will fit in a fixnum,
  763. ;;; so all single-floats larger than most-positive-fixnum can be precisely
  764. ;;; represented by an integer.]
  765. ;;;
  766. (defun %unary-round (number)
  767.   (number-dispatch ((number real))
  768.     ((integer) number)
  769.     ((ratio) (values (round (numerator number) (denominator number))))
  770.     (((foreach single-float double-float))
  771.      (if (< (float most-negative-fixnum number)
  772.         number
  773.         (float most-positive-fixnum number))
  774.      (truly-the fixnum (%unary-round number))
  775.      (multiple-value-bind (bits exp)
  776.                   (integer-decode-float number)
  777.        (let* ((shifted (ash bits exp))
  778.           (rounded (if (and (minusp exp)
  779.                     (oddp shifted)
  780.                     (eql (logand bits
  781.                          (lognot (ash -1 (- exp))))
  782.                      (ash 1 (- -1 exp))))
  783.                    (1+ shifted)
  784.                    shifted)))
  785.          (if (minusp number)
  786.          (- rounded)
  787.          rounded)))))))
  788.  
  789.  
  790. (defun rational (x)
  791.   "RATIONAL produces a rational number for any real numeric argument.  This is
  792.   more efficient than RATIONALIZE, but it assumes that floating-point is
  793.   completely accurate, giving a result that isn't as pretty."
  794.   (number-dispatch ((x real))
  795.     (((foreach single-float double-float))
  796.      (multiple-value-bind (bits exp)
  797.               (integer-decode-float x)
  798.        (if (eql bits 0)
  799.        0
  800.        (let* ((int (if (minusp x) (- bits) bits))
  801.           (digits (float-digits x))
  802.           (ex (+ exp digits)))
  803.          (if (minusp ex)
  804.          (integer-/-integer int (ash 1 (+ digits (- ex))))
  805.          (integer-/-integer (ash int ex) (ash 1 digits)))))))
  806.     ((rational) x)))
  807.  
  808.  
  809. (defun rationalize (x)
  810.   "Converts any REAL to a RATIONAL.  Floats are converted to a simple rational
  811.   representation exploiting the assumption that floats are only accurate to
  812.   their precision.  RATIONALIZE (and also RATIONAL) preserve the invariant:
  813.       (= x (float (rationalize x) x))"
  814.   (number-dispatch ((x real))
  815.     (((foreach single-float double-float))
  816.      ;; Thanks to Kim Fateman, who stole this function rationalize-float
  817.      ;; from macsyma's rational. Macsyma'a rationalize was written
  818.      ;; by the legendary Gosper (rwg). Gosper is now working for Symbolics.
  819.      ;; Guy Steele said about Gosper, "He has been called the
  820.      ;; only living 17th century mathematician and is also the best
  821.      ;; pdp-10 hacker I know." So, if you can understand or debug this
  822.      ;; code you win big.
  823.      (cond ((minusp x) (- (rationalize (- x))))
  824.        ((zerop x) 0)
  825.        (t
  826.         (let ((eps (if (typep x 'single-float)
  827.                single-float-epsilon
  828.                double-float-epsilon))
  829.           (y ())
  830.           (a ()))
  831.           (do ((xx x (setq y (/ (float 1.0 x) (- xx (float a x)))))
  832.            (num (setq a (truncate x))
  833.             (+ (* (setq a (truncate y)) num) onum))
  834.            (den 1 (+ (* a den) oden))
  835.            (onum 1 num)
  836.            (oden 0 den))
  837.           ((and (not (zerop den))
  838.             (not (> (abs (/ (- x (/ (float num x)
  839.                         (float den x)))
  840.                     x))
  841.                 eps)))
  842.            (integer-/-integer num den))
  843.         (declare ((dispatch-type x) xx)))))))
  844.     ((rational) x)))
  845.